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ABSTRACT 

We present A^-body simulations of planetary accretion beginning with 1 km radius planetesimals in orbit about 
a 1 Mq star at 0.4 AU. The initial disk of planetesimals contains too many bodies for any current A^-body code 
to integrate; therefore, we model a sample patch of the disk. Although this greatly reduces the number of bodies, 
we still track in excess of lO*" particles. We consider three initial velocity distributions and monitor the growth of 
the planetesimals. The masses of some particles increase by more than a factor of 100. Additionally, the escape 
speed of the largest particle grows considerably faster than the velocity dispersion of the particles, suggesting 
impending runaway growth, although no particle grows large enough to detach itself from the power law size- 
frequency distribution. These results are in general agreement with previous statistical and analytical results. We 
compute rotation rates by assuming conservation of angular momentum around the center of mass at impact and 
that merged planetesimals relax to spherical shapes. At the end of our simulations, the majority of bodies that have 
undergone at least one merger are rotating faster than the breakup frequency. This implies that the assumption of 
completely inelastic collisions (perfect accretion), which is made in most simulations of planetary growth at sizes 
1 km and above, is inappropriate. Our simulations reveal that, subsequent to the number of particles in the patch 
having been decreased by mergers to half its initial value, the presence of larger bodies in neighboring regions of 
the disk may limit the validity of simulations employing the patch approximation. 

Subject headings: 



1. INTRODUCTION 

The "planetesimal hypothesis" states that planets grow 
within circumstellar disks via pairwise accretion of small solid 
bodies known as planetesimals (Chamberlin 1905; Safronov 
1969; Hayashi et al. 1977). The process of planetary growth 
is generally divided for convenience and tractability into sev- 
eral distinct stages. In the first stage, microscopic grains col- 
lide and grow via pairwise collisions while settling towards the 
midplane of the disk. If the disk is laminar, then the solids 
may collapse into a layer that is thin enough for gravitational 
instabilities to occur (Edgeworth 1949; Safronov 1960; Gol- 
dreich and Ward 1973; Youdin and Shu 2002; Garaud and Lin 
2004; Johansen et al. 2007); such instabilities would have pro- 
duced planetesimals of ^ 1 km radius at 1 AU from the Sun. If 
the disk is turbulent, then gravitational instabilities may be sup- 
pressed because the dusty layer remains too thick. Under such 
circumstances, continued growth via binary agglomeration de- 
pends upon (currently unknown) sticking and disruption prob- 
abilities for collisions between larger grains (Weidenschilling 
and Cuzzi 1993; Weidenschilling 1995). The gaseous compo- 
nent of the protoplanetary disk plays an important role in this 
stage of planetary growth (Adachi et al. 1976; Weidenschilling 
1977; Rafikov 2004). 

Once solid bodies reach kilometer-size (in the case of the 
terrestrial region of the proto-solar disk), gravitational interac- 
tions between pairs of solid planetesimals provide the domi- 
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nant perturbation of their basic Keplerian orbits. Electromag- 
netic forces, collective gravitational effects, and in most cir- 
cumstances gas drag, play minor roles. These planetesimals 
continue to agglomerate via pairwise mergers. The rate of solid 
body accretion by a planetesimal or planetary embryo is deter- 
mined by the size and mass of the planetesimal/planetary em- 
bryo, the surface density of planetesimals, and the distribution 
of planetesimal velocities relative to the accreting body. The 
evolution of the planetesimal size distribution is determined by 
the gravitationally enhanced collision cross-section, which fa- 
vors collisions between bodies having smaller relative speeds. 
Runaway growth of the largest planetesimal in each accretion 
zone appears to be a likely outcome. The subsequent accumu- 
lation of the resulting planetary embryos leads to a large de- 
gree of radial mixing in the terrestrial planet region, with giant 
impacts probable. Growth via binary collisions proceeds until 
the protoplanets become dynamically isolated from one another 
(Lissauer 1987, 1995). 

Numerous groups have attempted to examine the accumula- 
tion and dynamics of 1 km planetesimals via numerical simula- 
tions. The statistical approach (Greenberg ef aZ. 1978; Kolvoord 
and Greenberg 1992) came to be known as the particle-in-a-box 
(PIAB) method. More recent PIAB investigations, also begin- 
ning with 1 km planetesimals, have been performed by Wei- 
denschilling et flZ. (1997) and Kenyon and Bromeley (2006). 
PIAB assumes that the velocity distribution of planetesimals 
is a smooth function. Each particle effectively sees a "sea" of 
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particles. The advantage of this model is that one need not sum 
the gravity between all the bodies; rather, the relative veloc- 
ities and impact parameters of two colUding particles can be 
randomly selected from a simple function. 

Lecar and Aarseth (1986) were the first to apply A'-body 
modeling of planetary growth with a simulation of 100 lunar- 
sized bodies that interacted directly and accreted. Com- 
puter power and algorithm sophistication increased through the 
1990's, but these improvements were not sufficient to enable 
the direct A^-body modeling of growth from 1 km planetesi- 
mals. Nevertheless, great strides were made toward understand- 
ing the final stages of planet formation {e.g., Agnor et al. 1999; 
Chambers 2001; Kokubo and Ida 2002; Raymond et al. 2006; 
O'Brien etal. 2006; Morishimaefa/. 2008). The highest (the 
number of particles) simulation to date is that of Richardson et 
al. (2000), which modeled one milhon 150 km radius particles 
for lO-' years. 

In this investigation we simulate 1 km planetesimal growth 
with an Af-body model. Our approach is not to consider an en- 
tire disk of planetesimals, as calculating gravitational interac- 
tions between more than 1 trilUon particles is an intractable task 
for the foreseeable future. Instead we focus on small, square, 
shearing patches of the disk, containing up to 10^ particles. In 
this way, we follow growth over several orders of magnitude, 
compare self-consistent and statistical calculations, and lay the 
groundwork for future investigations of direct simulations of 
planetesimal growth. 

In §2, we describe the numerical integration techniques. In 
§3, we summarize the initial conditions of our simulations. In 
§4, we present the results for our baseline model, in which the 
magnitude of the initial velocity dispersion is equal to the es- 
cape speed of 1 km planetesimals. In §5, we compare the results 
of two simulations with different initial velocity dispersions. In 
§6, we discuss some of the key results to emerge from these 
simulations. Finally, in §7, we draw more general conclusions, 
extrapolate our results to longer times and larger orbital radii, 
and describe future directions of research. Appendix A lists all 
symbols and abbreviations used in this article. Appendix B re- 
views an analytic method for approximating planetesimal accu- 
mulation. Appendix C presents three of our simulations to 2000 
orbits; these results have been relegated to an appendix because 
they probably suffer from systematic errors due to the small 
physical size of the region being simulated combined with the 
flat slope of the size-frequency distribution of planetesimals at 
this epoch. 

2. NUMERICAL TECHNIQUES 

We use flie code PKDGRAV (Richardson et al. 2000; Stadel 
2001 ; Wadsley et al. 2004) to perform the integrations. This is a 
parallel, highly scalable A'-body algorithm, originally designed 
for cosmological simulations (see, e.g., Moore et al. 1998). The 
code incorporates features such as multipole expansions, mul- 
tistepping, and binary merging to increase speed (described be- 
low). Despite these sophisticated techniques, PKDGRAV can- 
not integrate more than ^10^ colliding 1 km planetesimals for 
10^ years in a reasonable amount of CPU time with our model. 
Since the collision timescale is of order minutes (see §2.3), 
our "baseline model", presented in §4, required over 4x10^ 
timesteps. We describe the equations of motion of the patch 
in §2.1, our collisional model in §2.2, the basic principles of 
PKDGRAV in §2.3, the numerical approximations necessary to 
complete these simulations in §2.4, and our methods of verifi- 



cation in §2.5. 



2.1. Equations of Motion 



Although carrying out the simulations in a patch greatly re- 
duces A', it adds the new complication of simulating Keplerian 
shear: we must account for the gradient in the Sun's gravita- 
tional field and the geometry of the disk. If we examine small 
patches of the disk such that W ^r, where W is the width and 
length of the square patches and r is the distance from the Sun, 
then we may approximate quantities, such as surface density, 
scale height, and circular velocity, as varying linearly in the 
patch. 

Two coordinate systems are used in this model: heliocentric 
cylindrical and locally Cartesian. The heliocentric system is the 
more physically meaningful system, as it is based on the geom- 
etry of the disk. In these coordinates, the origin is located at 
the position of the central star, and (r, 6, z) have their standard 
meanings. We implement the Cartesian system inside the patch, 
with the origin at the center of the patch. For small 9, these two 
coordinate systems are related by the foUowing expressions: 
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where r patch is the heliocentric radius of the center of the patch 
and Q.patch is the Keplerian orbital frequency of the patch. Es- 
sentially the radial direction is mimicked by x and the tangential 
direction by y. Similarly, the velocities are related by 
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The z-velocities and positions are identical in the two coordi- 
nate systems. 

For the purpose of the integration, we use the Cartesian sys- 
tem, as outlined in Wisdom and Tremaine (1988), which was 
based on Hill (1878). The equations of motion inside the patch 
are 

X- 2Q.patchy-^^patch^ = 

y+2flpatchX = (3) 

''''patch'" ~ dz' 

where <j> represents the interparticle potentials in the disk. Wis- 
dom and Tremaine assume massless particles (Vcj) = 0), and 
provide the general solution 

X = Xg+A COSiflpatchO+Bsiniilpatcht), 

y = yg- l^patchXgt - 2A &m(npa,cht) + 2Bcos{Qpatcht), (4) 

Z = CcOS(flpatcht) +Dsin(^lpatcht), 

where A,B,C, and D are arbitrary constants that correspond to 
the amplitude of the excursions from a circular, planar orbit, 
and Xg and yg are the guiding centers of the epicycles of the par- 
ticles. The jftpatchXg term is known as the "shear" and repre- 
sents the Keplerian motion relative to that at x = 0. For negative 
X, the shear is positive, as these particles are closer to the star. 
The y-motion is the tangential motion relative to the center of 
the patch, not the azimuthal motion of the patch as a whole. 

The equations of motion are integrated with a second-order 
generalized leapfrog scheme derived via the operator splitting 
formalism of Saha and Tremaine (1992). For the case of or- 
bits in the patch, this integrator prevents secular changes in the 
guiding center of the epicycle. 
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2.2. Collision Model 

Our model assumes perfect accretion. In our implementa- 
tion, collisions result in one perfectly spherical particle, and all 
particles have the same density. To conserve angular momen- 
tum, the resultant particle receives the net angular momentum 
of the two colliders (spin plus motion relative to the center-of- 
mass of the colliding particles), without dissipation (although 
energy is decreased by 20% during a collision). Our assumption 
of perfect accretion allows particles to spin faster than break-up, 
the latter being given by 



^Jll'K/Gppl. 



(5) 



(For a density of 3 g/cm', p„in = 1.9 h.) Real gravitational ag- 
gregates do not merge completely under such circumstances, 
and hence this assumption will tend to artificially increase plan- 
etesimal masses. Conversely, spheres have the smallest possible 
cross-section for a given volume, so our algorithm is suppress- 
ing growth with this approximation. Therefore these two fea- 
tures of our model counteract each other, although we cannot 
say to what degree they cancel each other out. 

The assumption that the particles are always gravitational 
aggregates may break down during high-energy collisions, in 
which case melting may occur. For tractability, we ignore this 
possibility. We also ignore any potential radiogenic heating, 
e.g., by ^^Al, as the melting timescale is much longer than the 
colUsional timescale. 

2.3. PKDGRAV 

PKDGRAV is very efficient because it calculates gravita- 
tional forces rapidly. The code works by recursively dividing 
the physical space, the "domain", into smaller cells that are or- 
ganized in a tree. The gravitational forces are then calculated 
by traversing this tree (Barnes and Hut 1986). To calculate the 
forces on a particle from a given cell, the code determines the 
apparent size of the cell, that is the angle the cell subtends as 
viewed from the particle. If this angle is smaller than some 
minimum angle, O, then all moments up to hexadecapole of 
the particles in the cell are used to calculate the force. Other- 
wise the cell is opened and the test is repeated on the subcells. 
This approximation changes the speed of the A^-body calcula- 
tion from OiN^) to OiNlogN). For further details see Stadel 
(2001) and Wadsley et al. (2004). For the simulations described 
here, we used O = 0.7 radians. 

PKDGRAV also permits the use of periodic boundary con- 
ditions, which are required to prevent the patch from self- 
coUapsing. To do this, each patch is reproduced 8 times, in 
the form of "ghost cells" that completely surround the actual 
patch in the x and y directions. Eq. (4), however, also requires 
that particles in the ghost cells move at the appropriate shear 
velocity. The centers of the ghost cells therefore move in the 
y-direction as prescribed by Eq. (4). 

Gravity calculations per time interval are minimized by 
a multistepping algorithm. Multistepping divides the base 
timestep into smaller intervals. For our simulations the timestep 
is based on the largest acceleration a particle feels at each 
timestep. In this algorithm isolated planetesimals move at the 
base timestep, those, set to about 100 steps per orbit. As particles 
approach each other their timesteps drop to a minimum value, 
until they either miss or collide. Each smaller timestep is a fac- 
tor of 2 shorter then the previous in order to keep all of the base 



gravitational kicks commensurate. Typically 95% of particles 
are on the longest timestep, but we still resolve all collisions. 

Two timescales are relevant in this problem, the orbital time 
and the crossing time of two particles that just miss. The or- 
bital period, P, is of order one year, but the crossing time is of 
order minutes. For two particles of radius R to just miss, on a 
parabolic orbit, the crossing time is 



■■ r/ve. 
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and is therefore independent of the masses. In this scenario and 
with densities of 3 g/crn^, the crossing time is 772 s. P and f„o.s.v 
determine those and tmin, the maximum and minimum possible 
timesteps. To be conservative (note that the version of PKD- 
GRAV employed in this investigation does not integrate Eq. [3] 
symplectically), we set these values as 



tbase = r]P (7) 



and 



tmin ~ V^crossj (8) 

where is a scale factor, chosen such that the integrals of mo- 
tion are constant to a satisfactory degree. Convergence tests 
showed that for the systems we consider, rj = I /300 provided 
the necessary accuracy. The number of available timesteps, (, 
is 

(9) 



C=l+log,(^). 
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At 0.4 AU and p = 3 g/cm', this translates to C = 14. 

2.4. Numerical Approximations 

Although PKDGRAV is fast relative to other integration 
techniques, we must make several additional approximations 
in order for the simulation to be practical, and to appropriately 
model Eq. (4). These approximations do not alter the physics 
appreciably, and resulted in nearly a factor of 10 speed-up. 

As mentioned above, the timestep for each particle is deter- 
mined by its current acceleration, which is a function of the 
local mass density. Particles on large timesteps are not close to 
other particles and cannot be close to collision. We therefore 
set a ceiling for timesteps to search for collisions. This small 
modification increases speed by a factor of 2 - 3 depending on 
the surface density of the patch (and hence the local density). 

Occasionally during the evolution of the patches, binaries 
(gravitationally bound pairs of particles) form. This is a nat- 
ural result of 3-body encounters. Although it is preferable to 
allow these binary systems to evolve normally, they often be- 
come stuck in the shortest timesteps, greatly diminishing the 
advantage of multistepping. Therefore, we artificially merge 
them. In our simulations, binary merging accounts for about 
0.01% of all merging events and therefore this small change in 
the total energy of the simulation should be negligible. Simula- 
tions of globular clusters have shown that tight binaries tend to 
get tighter during a 3-body encounter, and loose binaries tend to 
become looser (Heggie 1975). To make binary merging as re- 
alistic as possible we therefore also require the eccentricity of 
the particles involved to be less than 0.9, so that wide binaries 
may still be disrupted. 

2.5. Verification 

Given the complexity of this problem (non-inertial frame, 
large range of growth, and numerical approximations), we need 
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to quantify the accuracy of our methodology. The patch frame- 
work provides some inherent tests, but we must also understand 
the statistics of the problem; this patch is supposed to be a rep- 
resentative piece of a much larger annulus of material. At some 
point, the number of particles drops to a sufficiently small num- 
ber that the results cannot be trusted. In this subsection we de- 
scribe our methodology for verifying the results presented in 
§§4-5. 

In the shearing model, the particles are in a non-inertial 
frame. Therefore the integrals of motion normally associated 
with dynamics (momentum, energy, angular momentum), are 
not directly appUcable. In this formalism there are two con- 
stants of motion we use to check the validity of these simula- 
tions. In the nomenclature of Wisdom and Tremaine (1988), 
these conserved integrals are 



ately. We parameterize this effect as 



EN 



(10) 



patch 



where M patch is the total mass in the patch. These parameters 
essentially correspond to the center of mass velocity. For this 
situation, the center of mass should remain motionless; we need 
to verify that u and w remain much less than the random veloc- 
ities. Note that our definition includes the masses of the parti- 
cles, whereas the Wisdom and Tremaine (1988) definition did 
not, as they weighted all particles equally. 

Eq. (10) would be sufficient if we were integrating particles 
without growth, which we are not. The patch model at least 
requires that no particle is dominant. Therefore a zeroeth or- 
der requirement is that no particle reaches a mass equal to that 
of the sum of the remaining particles. This constraint, how- 
ever, fails to take into account how the most massive particles 
alter the dynamics of the disk. As a typical particle passes by 
the largest mass particle, it receives a kick, and energy associ- 
ated with Keplerian motion may be transformed into the ran- 
dom motions of the swarm, increasing the velocity dispersion, 
i.e., "viscous stirring" (Wetherill and Stewart 1989). Afterward 
that increase in random velocity should be damped down by 
subsequent interactions with other typical particles. To verify 
that no particles grow so large as to dominate the stirring in the 
patch, we consider the ratio of the gravitational stirring of the 
largest particle to the sum of all other particles: 



w 



(12) 



where W is the patch width. The equations of motion do not de- 
pend on eccentricity, so there are no numerical problems if this 
occurs, only concerns about the physical interpretation of this 
phenomenon. Our choice of W is very small compared to the 
size of the terrestrial annulus. As long as (3 remains less than 
or close to unity, the radial mixing throughout the disk is neg- 
ligible, and our interpretations are independent of the particles' 
eccentricity. Only if (3 grows to large values is there cause for 
alarm. We focus on the second largest /3 value, as occasional 
strong kicks could temporarily result in large eccentricities for 
one particle. 

3. INmAL CONDITIONS 

We assume a surface density of 37.2 g/cm^ at 0.4 AU, 30% 
higher than predicted by the minimum mass solar nebula model 
(Hayashi 1981), but consistent with previous studies (Richard- 
son et al. 2000). The bulk density of the planetesimals, Ppi, is 
taken to be 3 g/cm^. The mass of a planetesimal is thus 
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Assuming the surface density of the disk scales as S oc a"', the 
number of planetesimals in the annulus 0.4 AU < a < 2.5 AU 
is 



rripi 



3.5 X 10 
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As we are unable to directly integrate trillions of particles, we 
must reduce A' to a tractable value by dividing the disk into 
patches. 

Most of our simulations begin with the root mean squared 
(RMS) velocity dispersion, Vrms, set to equal to the escape 
speed of 1 km planetesimals, Vesc (Safronov 1969; Stewart and 
Wetherill 1988). This relationship assumes that the system 
is relaxed, which may or may not be the case depending on 
the processes and timescale to form planetesimals. This RMS 
speed is the magnitude of the random motion of the particles, 
but it is not distributed isotropically due to the 2:1 axis ratio of 
the epicycle; instead the distributions are described by 



E 



/<;,„ 



Nim]-ml^ 



(11) 



(Lissauer & Stewart 1993). Exactly how large S can grow is not 
clear a priori, but values below ~ 0.1 should be satisfactory. 

Should the mass distribution in the patch suggest that sig- 
nificantly massive planetesimals are present in the disk, but not 
present in a patch, then the patch is not modeling a large enough 
region of the disk. We check this possibility by plotting rn^Nk as 
a function of m, where Nk is the number of particles in mass bin 
k = mjm\. We will see in §4.4.2 and Appendix C that unmod- 
eled large bodies are hkely to become a problem as our model 
evolves. 

A final point of concern is the size of the epicycles of parti- 
cles compared to the size of the patch. Should the eccentricity 
of a particle grow large enough that the radial excursions (2ae, 
where a is semi-major axis, e is eccentricity) exceed the size 
of the patch, then we may not be sampling the region appropri- 



3 ^esc 



(15) 



(Binney and Tremaine 1994). For a radius of 1 km and ppi = 3 
g/cm^ particles, Vesc is 1.29 m/s. 

From this velocity distribution, we calculate the equilibrium 
vertical density profile of the disk. To do this we assume the 
disk is in a state of hydrostatic equihbrium; the "pressure" of 
the vertical component of the velocity distribution maintains its 
thickness. The vertical component of solar gravity increases 
linearly with distance from the midplane, therefore the density 
follows a Gaussian profile, 



P = Poe 



(16) 



where p is the density, po is the density at the midplane, and Zq 
is the scale height of the disk. The (Gaussian) scale height is 
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determined by the z-velocity distribution, the mass of the cen- 
tral star and the orbital radius. Hydrostatic equiUbrium allows 
a calculation of the scale height: 




(17) 



The scale height also depends on the self-gravity of the plan- 
etesimal disk. To compensate for this, we implement a "vertical 
frequency enhancement" to simulate an infinite plane sheet of 
matter so that 



n, = n 
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(18) 



The second term in Eq. (18) is an analytic restoring force that 
simulates the disk's self-gravity and reduces Zq. For our base- 
line model (Simulation Li presented in §4), Cl^ is 0.2% larger 
than ^patch- 

Our nomenclature for simulations is based on the relative 
size and the initial dynamical state. L stands for "large", M for 
"medium" and S for "small". The subscript is the ratio of the 
initial velocity dispersion to the escape speed of a 1 km plan- 
etesimal. In Table 1, we present the initial conditions for these 
simulations. 

We assume that scattering is very effective, and that gas drag 
is negligible. Simulation Li, Mi and Si presume the planetesi- 
mals are in a form of equilibrium. However, there remain many 
unknowns in the formation of these planetesimals and varia- 
tions of a factor of 2 are plausible. We therefore ran two inte- 
grations with non-equilibrium initial conditions. We integrated 
one system that began with vrms = 0.5vesc (Simulation Mq.s), 
and one with vrms = '^^esc (Simulation M2). In these simulations 
the scale height was set by Eq. (17). These two simulations are 
presented in §5. 

Each of the simulations has been run at least until the num- 
ber of bodies in the patch has been halved, fi/2. This time is 
related to the mean free time, r, between coIUsions, which can 
be approximated as 

T^^^ ~?l/2, (19) 

where a pi is the gravitational cross-section, o'(l + Vgj^./v^_5), 
and po is given by 

P^ = ^:r- (20) 
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Eq. (19) is only valid at the midplane at time f = 0. The density 
changes with z, and vrms change with time. Nonetheless, this 
simple model provides valuable insights into the dynamics. For 
the initial conditions we use, r k, 500 orbits. 

Richardson et al. (2000) modeled a planetesimal disk with 
10^ particles, but they simulated a later stage of planetary 
growth in which the collision rate was much lower than it is 
here. Therefore, our integrations proceed slower, even with the 
advanced methodologies described §2.3. If we set Npatch 10^, 
then the simulations are large, but we may still examine sev- 
eral different initial conditions. We choose a square patch of 
width W = 10"^ 

f patch for our largest simulation, Li. In our 
model, this choice corresponds to an initial A'^ of 106,130, a 
small enough number to be tractable, but large enough to be 
statistically significant. Note, however, that the medium-sized 
simulations have dimensions W = 5 x 10~'^r, and the small one 
W = 2.5x I0~^r. These cases run faster, and are included to test 
various assumptions in the baseline case. 



4. THE BASELINE MODEL 

The simulations of our baseline model (those with subscript 
1) begin with 1 km radius particles with an RMS velocity equal 
to that of the escape speed of 1 km particles with a density of 
3 g/cm^. From this dispersion the initial scale height is set by 
Eq. (17). The simulation with the most particles, Li, required 
354 orbits (89.6 years) to halve the total number of particles. 
The Ml and Si simulations required 349 and 361 orbits, re- 
spectively, to reach ti/2- In this section we describe the results 
of these three simulations. Table 2 summarizes some of the re- 
sults of our simulations at ti/2. At this time, fi/2, the patch is 
nearing a condition in which it is too small to correctly model 
the velocity distribution (see §4.4.2); results of these simula- 
tions beyond ti/2 are presented in Appendix C. 



4.1. Mass Growth 

The mass spectra of Simulations Li, Mi, and Si at 100, 200, 
300 and 354 orbits (fi/2 for Li) are plotted in Fig. 1. Figure 2 
shows the (differential) mass distributions with logarithmic bins 
of size 2'oti,/ = 1,2,3, i.e., the first bin contains all particles 
of mass mi, the second bin contains particles of 2mi and 3mi, 
the third bin contains particles in the range 4m 1 < m < 7m 1, 
etc. The histograms shown in Fig. 2 appear to follow a power 
law size distribution. The one apparent exception is the sin- 
gle particle in the largest bin for Li at ?i/2, which is detached 
by one bin from the remaining particles (the "swarm"). How- 
ever, upon closer inspection, this particle has a mass of 275 mi, 
whereas the second and third largest particles have masses of 
122 and 116 mi, respectively. All three of these particles just 
miss being members of the 128 - 255 mi bin. Therefore, the 
distributions presented in Fig. 2 suggests that runaway growth 
has not occurred. 




10 100 

Mass (mj 



10 100 
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Fig. 1 . — The semi-log (differential) mass distributions for Simula- 
tions Li (solid line), Mi (dashed line) and Si (dotted line) after 100 
orbits (top left), 200 orbits (top right), 300 orbits (bottom left) and 354 
orbits (bottom right). 
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Fig. 2. — The log-log (differential) mass distributions of Simulations 
Li (solid line), Mi (dashed line) and Si (dotted line) at 100 orbits (top 
left), 200 orbits (top right), 300 orbits (bottom left), and 354 orbits 
(bottom right). Note that the Mi simulation data have been offset to 
the right by 0.1 and the Si data by 0.2 in order to improve readabiUty. 

To further examine the nature of the mass distribution, we 
matched the mass distribution to single parameter fits. We con- 
sider a power law, 



N{k) oc r 



and an exponential, 



Nik) oc e 



-k/c 



(21) 



(22) 



We fit the constants b and c by using the line method to min- 
imize (Press et al. 1996) with one sigma uncertainties as- 
sumed to be \/Nk + .75 + 1 . For small Nk this gives a better ap- 
proximation to the distribution of for sparsely sampled data 
(Gehrels 1986), e.g., some bins at large k have zero particles. 
The fits are constrained so that the total mass is the same as 
in the simulation, and the number of degrees of freedom is ap- 
proximately equal to m,nax/nii . We must take into consideration 
the fact that the first bin is unusual in that all the particles were 
originally in that bin. Therefore we also consider fits that ex- 
clude the data at m = m\. These fit parameters are denoted with 
a prime. We found that exponential fits are always poor, so we 
will focus on the power law fits. 

In Fig. 3 we show the mass distributions of Simulations Li 
(top). Ml (middle) and Si (bottom), with the analytic fits for 
comparison at each simulation's ti/2- The values of all the pa- 
rameters are presented in Table 3. In that table we see that the 
values of unreduced are all significantly better when the first 
data point is excluded. All fits include the large-mass particles 
and we conclude that these largest particles do not represent a 
new class (i.e., runaway growth has not occurred). Note that 
forcing the fitted distribution to contain the same total mass as 
the simulation forces the exponential cases to converge to poor 
fits. 



Fig. 3. — Comparison of the mass distributions at ti/2 for the baseline 
simulations. The histograms represent the mass distribution at each 
model's value of ti/2 and the straight lines are the power law fits, soUd 
includes the first bin, dashed does not. 

In Fig. 4 we compare the merger rate (-dN/dt) of Simula- 
tion Li with that predicted by the mean free time of the patch. 
The prediction comes from the following relation 



dN _ Npo < a pi > vrms 
dt <m> ' 



(23) 



where the brackets indicate quantities averaged over the entire 
swarm. We assume the cross-section to be the gravitationally 
enhanced area of two particles with the average radius. The 
right hand side of Eq. (23) is the number of particles divided 
by the mean free time. Initially the mean free time argument 
underpredicts the merger rate and later it is about right. Over- 
all this simplified estimate of the merger rate remains near the 
actual value. The rate is dropping because the total number of 
particles in the patch is decreasing due to merging. 
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Fig. 4. — The evolution of tiie merger rate of Simulation Li. The 
solid line is that observed in the patch (sampled once per orbit), and 
the dashed line is that predicted by the mean free time of the patch (re- 
calculated each orbit from the instantaneous values of vrms, < Upi >, 
po and < m >) in Simulation Li. 

Next we compare our results to that of the constant, linear, 
and product solutions to the coagulation equation (Wetherill 
1990; see also App. B) in Simulation Li. These three solu- 
tions are relatively simple compared to more recent derivations 
(see e.g., Kenyon and Luu 1998; Kenyon and Bromley 2004). 
However, the recent models do not provide analytic solutions 
for the number of particles in each mass bin (see Eqs. [B9], 
[BIO] and [B12]). We find the three collisional probability co- 
efficients (see Eqs. [B6 - B8]) have values of i^i = 5.3 x 10"**, 
1^2 = 1.84 X 10"^ and I'i = 2.65 x 10"^ for Simulation Li. 

In Fig. 5 we plot the growth of some of the largest masses 
in Simulation Li through 354 orbits. Also shown are the pre- 
dicted largest masses from the three solutions to the coagula- 
tion equation, with collisional probabilities that are constant as 
a function of time, and equal to the values listed above. The 
growth in our A^-body model follows the product solution to 
the coagulation equation for about 250 orbits, but then the two 
diverge as the A^-body model predicts faster growth. This di- 
vergence is probably a result of the product solution's failure to 
conserve mass for f > ti/2- At fi/2 the largest particle has a mass 
of 276 nil ■ Nonetheless, the agreement over the majority of the 
simulation suggests that even the Wetherill (1990) product co- 
agulation model is a reasonable representation of early growth 
of planetesimals. 
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Fig. 5. — The evolution of the largest, second largest, fifth largest and 
tenth largest masses in Simulation Li. The solid red line corresponds 
to nimax in the product solution to the coagulation equation, green to 
the linear, and blue to the constant. 

In Fig. 6 we further compare Simulation Li to the Wether- 
ill (1990) product model. The six panels examine six different 
ranges of k. For mass ranges, the predicted number of particles 
is determined by summing all mass bins in that range. From 
this figure we see that no solution to the coagulation fits the data 
over all values of k. At low masses (k < 10), the linear solution 
is the best-fit to the data, but at larger masses the product solu- 



tion is better. However, at these larger values of k, the product 
solution still differs from the actual distribution by more than a 
factor of 2. 
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Fig. 6. — A comparison of the distribution of masses in Simulation 
Li with those predicted by coagulation theory for 6 different mass 
ranges: k = 1 (bottom left), yt = 10 (middle left) 20 < < 50 (mid- 
dle left), 50 < *r < 100 (middle right), 100 < A: < 200 (top left) and 
200 <k< 300 (top right). 



4.2. Velocity Dispersion 

Figure 7 shows the evolution of the RMS velocity of all par- 
ticles, the escape speed of the largest particle, and the escape 
speed of the average-mass particle in Simulation L]. The ve- 
locity dispersion slowly grows to a value of 2 m s"' at fi/2 
(which is evidence of viscous stirring), while the escape speed 
of the largest particle grows to a value of 8.4 m s"'. At fi/2 
the largest body has a mass of m,„a^ = 276m i and a radius of 
6.5 km. From these values we can identify the gravitational fo- 
cusing factor for the largest body, the ratio of its gravitational 
cross-section to its geometrical cross-section. This factor is 



Fg = l + 



2Gm„ 



R 



= l+(^) 



Vrms 



(24) 
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Fig. 7. — The evolution of the RMS velocity of the patch (solid line), 
the escape speed of the largest particle (dashed line) and the escape 
speed of the mean particle (dotted line) in Simulation Li. 
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Fig. 8. — Evolution of the gravitational focusing factor of the most 
massive particle in Simulations Li (sohd line), Mi (dashed hne), and 
Si (dotted line). In larger simulations the focusing grows larger be- 
cause the biggest particles are more massive, but the velocity disper- 
sions are about equal in the three simulations (see Table 2). 



We evaluate the ratios of crms to sin irms in Fig. 9. (Note the 
inclinations are in a regime such that the difference between sin 
; and i is about 1 part in 10^".) The ratios begin at 2.35 but drop 
to ~ 2. 15, similar to previous results (Greenzweig and Lissauer 
1990; Kokubo and Ida 1996). 
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Fig. 9. — Evolution of eRus/ sini^s for Simulations Li (solid Hne), 
Ml (dashed line), and Si (dotted line). The values drop from ^2.35 
to ^2.15 over ~200 orbits. At that time the ratios level out, which 
suggests the initial conditions were not in equilibrium. 

We investigate the role of dynamical friction by examining 
cgMS and Irms as a function of mass at /i/2- For each popu- 
lated bin, we computed these values (if only one body occu- 
pied a bin we used its e and ; values), and show the results 
in Fig. 10. As expected, larger masses are dynamically colder 
than smaller masses, however they have greater energy, imply- 
ing that equipartition of energy is not achieved. 




Mass (mj 

Fig. 10. — The values of brms and Irms as a function of mass at ti /2 for 
the three equilibrium simulations. Filled squares represent eccentric- 
ity, open triangles inchnation. In the Li simulation, the largest mass 
planetesimal (275mi) is not shown; its values are e = 2.6 x 10"* and 
i = 1.8 X 10~'. For reference, the dashed line represents equipartition 
of energy in e, dotted in /, normalized to the values for k = 10. 



4.3. Rotation Rates 
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We plot the rotation periods of all planetesimals in Simula- 
tion Li with m > mi at f in Fig. 11. Only fourteen bodies 
had periods longer than 1 day. Recall that initially all planetesi- 
mals have no spin. We plot the spin distribution over three mass 
ranges: m = 2m\, 3oti <m< lOmi, and m > lOmi. The peak 
period is at 1.05 hours, which is less than the minimum period 
for a spherical gravitational aggregate of density 3 g/cm^, see 
Eq. (5). Thus, this distribution suggests that our assumption 
(which is the standard one) of completely inelastic collisions 
resulting in mergers overestimates planetesimal growth rates. 
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Fig. 1 1 . — The final spin (at fi/2) periods for the planetesimals in Sim- 
ulation Li . The dashed vertical line represents the minimum period of 
a gravitational aggregate. 
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Fig. 12. — The evolution of the constants of motion for the nominal 
patch for Simulation Li . The values of u and w vary at a level 3 orders 
of magnitude below that of the random motions. 

4.4.2. Statistical Accuracy 

Although Simulation L| contains a large number of parti- 
cles (relative to modern A^-body simulations), it still represents 
a very small fraction of the terrestrial annulus. We therefore 
must characterize the robustness of our results. The most criti- 
cal aspect of this experiment is the mass of the largest particle. 
Should a particle reach a large enough mass that it inappropri- 
ately dominates the dynamics of the patch, then our assump- 
tions have broken down. 

First we consider /?, Eq. (12), which measures the radial ex- 
cursions of particles, see §2.5. We find that after 354 orbits of 
Simulation Li that (3max-\ (the second largest value) is still well 
below unity (Fig. 13). Therefore, the patch size for Simulation 
Li passes this requirement. 



0.35 



4.4. Accuracy Tests 

In this subsection we quantify the accuracy of our results (see 
§2.5). §4.4.1 measures the numerical accuracy of our code, and 
§4.4.2 describes the statistical accuracy of our method. 



4.4.1. Numerical Accuracy 

In Fig. 12 we examine the validity of Simulation Li by plot- 
ting u and w. These parameters measure the center-of-mass 
motion of the patch. Both u and w remain less 0.1 cm s~\ 
about 1000 times smaller than the typical random velocities in 
the patch and about 10^ times smaller than the shear across the 
patch, = 4720 cm s"' . We therefore conclude that this vari- 
ation is tolerable (Wisdom and Tremaine 1988). 
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Fig. 13. — The growth of the five largest values of (3 in Simulation Li 
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as a function of time. All values lie well below unity, indicating that 
no particle's epicycle is larger than the width of the patch. 



In Fig. 14, we plot the evolution of the stirring efficiency S 
(Eq. [11]) in the three baseline simulations. At fi/2 the largest 
particle is about one-eighth as effective at stirring as the rest 
of the swarm in Simulation Li, about one- seventh in Mi, and 
nearly one- fifth in Si. These values suggest we are nearing 
a situation in which the statistical accuracy of this simulation 
cannot be confirmed, but such a situation has not occurred yet. 
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Fig. 15. — The stirring power of particles as a function of mass at 100 
orbits (top left), 200 orbits (top right), 300 orbits (bottom left), and 354 
orbits (?i/2 for Simulation Li; bottom right). The soUd Une represents 
Simulation Li, dashed Mi, and dotted Si. After 354 orbits, the slopes 
are roughly flat and shows that our patches might not be large enough 
to contain a statistically significant number of particles. Note that the 
Simulation Mi data are offset by 0.1, and the Si data by 0.2. 



Fig. 14. — The stirring efficiency of the largest particles in Simula- 
tion Li (soUd line), Mi (dashed line), and Si (dotted line) relative to 
the rest of their respective swarms. 



Next we evaluate the distribution of rri^Nk as a function of 
mass. This quantity measures how effective each mass bin is at 
stirring the patch. As described in §2.5, the distribution can re- 
veal the possibility of large-mass bodies beyond the boundaries 
of the patch that could significantly change the dynamical char- 
acter of the planetesimal swarm. For this model, after about 
fi/2 there is a significant chance that large perturbers are close 
enough to affect the patches' dynamics. In Fig. 15 we plot m^Nk 
vs. log m at four times. If rrp-Nk is decreasing as a function 
of mass near the upper end of the mass distribution, then the 
patch is large enough to be a statistically representative piece 
of the annulus, and the velocity distribution of the patch should 
be close to that of the disk. If, alternatively, the distribution 
is increasing, then there exist nearby large particles that could 
dominate the stirring. We see at fi/2 the distribution is flat, sug- 
gesting we have reached the limit of our model. These curves 
suggest that results at later times may be inaccurate, hence our 
decision to exile those data in Appendix C. 



5. ALTERNATIVE VELOCITY DISTRIBUTIONS 

In this section we describe the results of two simulations that 
begin with different initial velocity dispersions than the base- 
Une model. The results of these simulations may be impor- 
tant since the initial velocity dispersion of planetesimals is ill- 
constrained. Earlier growth may occur too rapidly for the ve- 
locity dispersion to equilibrate (vrms = Vesc)- We thus explore 
a range of initial velocity dispersions. These simulations also 
show how sensitive the results of §4 are to variations in the ini- 
tial velocity dispersions. In Simulation M0.5 the initial velocity 
dispersion is set to 0.5v„x, in Simulation M2 it is set to Ivesc- 

For Simulations M0.5 and M2, fi/2 occurred after 253 and 542 
orbits, respectively. Figures 16 and 17 show the mass distribu- 
tions at fi/2 for each of these models. As expected, when the 
velocity dispersion is smaller, accretion proceeds faster, due to 
the increased gravitational focusing. Moreover, the largest par- 
ticles have a greater accretion advantage and so are more mas- 
sive at ?i/2 in M0.5 than in M2. 
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Fig. 16. — The mass functions of the non-equihbrium runs as a func- 
tion of time. The earliest time is offset to the right by 0.05, the next 
time is offset by 0.1, etc. Larger initial velocity dispersion suppresses 
runaway growth. These plots can be compared to Fig. 2. 

In Fig. 17 we show how our power-law fits (both with and 
without the A: = 1 data) to the log-log mass distribution, Eq. 
(21), compare to the actual distributions at each Simulation's 
fi/2- As with the baseline models (see Fig. 3), the exponential 
has a much larger unreduced value. The best fit parameters 
for these two simulations at time ?i/2 are Usted in Table 3. 
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Fig. 17. — Comparisons of the mass distributions at fi/2 for the Sim- 
ulations M2 (top) and M0.5 (bottom) and the best power law fits to the 
data, Eq. (21). The histograms are the mass distributions from the sim- 
ulations, the straight lines are the fits (solid includes first bin, dashed 
does not). 

We continue by plotting the evolution of the largest mass par- 
ticle for each non-equilibrium run in Fig. 18. As in Fig. 5, we 
also include the predictions of the Wetherill (1990) model. The 
linear solution is a good fit to the distribution of Simulation 
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M2, while the product solution appears to be a good fit to that 
of Simulation M0.5. 
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Fig. 18. — The growth of the largest mass in each model as a func- 
tion of time, and the predictions of the constant, linear, and product 
solution of the coagulation equation. 

In Fig. 19 we plot the evolution of the velocity dispersions. In 
the top panel (Simulation M2), we see that the velocity disper- 
sion actually drops initially, but then increases similarly to the 
other runs. This decrease results from the inelastic nature of the 
collisions, as well as, to a lesser degree, the deposition of trans- 
lational kinetic energy into rotational kinetic energy. However, 
we interpret this result with caution due to potential inconsis- 
tencies in the perfect accretion model (see §2.2). When the ini- 
tial velocity dispersion is too low, the dispersion quickly rises 
(compare to Fig. 7). The final values of Fg are 17.8 and 4.8 
for Simulations M0.5 and M2, respectively. As with the base- 
Une models, equipartition of energy has not occurred in these 
simulations. 
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equilibrium simulations. In Simulation M2 (top panel), the dispersion 
actually drops toward equilibrium. The final values of vrms are close 
to 2 m s"' , as in the baseline case. 

In Fig. 20 we plot the final spin distributions for the non- 
equilibrium patches. The peaks lie below the minimum gravi- 
tational aggregate period. Simulation M2 consists of especially 
fast rotators, due to a larger amount of kinetic energy of random 
motion available for transformation into rotational energy. 
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Fig. 20. — The final spin distributions for merged particles for the 
non-equilibrium simulations. The vertical dashed line is the minimimi 
period as defined by Eq. (5). For Simulation M2 very few particles 
have periods in excess of 10 hours, but Simulation M0.5 has several 
long-period planetesimals. 

6. DISCUSSION 

Several results stand out in §§4-5. First is that runaway 
growth has not begun to occur for any particle in any simula- 
tion. Second, the velocity dispersions of the patches remain 
close to the escape speed of the average-mass particle. Third, 
the growth rate is moderately sensitive to the initial velocity 
dispersion; changes of a factor of 2 in the initial RMS velocity 
can excite or retard the early growth rate of the most massive 
particles. Fourth, our collision model (§2.2) is too simplified in 
that it assumes spherical particles experiencing perfect accre- 
tion. Fifth, some aspects of the coagulation equation model the 
growth well, but discrepancy of at least a factor of 2 are present, 
and the best model (linear or product) appears to depend on the 
initial value of vrms- Sixth, the power law fits represent a real- 
istic model of the actual mass distribution. 

We present a summary of the results in Tables 2-3. Compar- 
ing Simulation Mi to the altemate-RMS velocity trial s, we see 
that its final properties lie between those of Simulations M0.5 
and M2. Therefore, m„u,x and fi/2 depend upon with the initial 
velocity dispersion in a systematic manner r. 

Growth proceeds easily in all our models, despite no initial 
seed. There is some indication from the final mass distributions 
that the annulus will develop particles with a mass in excess 
of lO'^mi by /i/2- By considering the stirring effects as a func- 
tion of mass we have found that by ?i/2 such large bodies could 
significantly modify the dynamical properties of our patches. 



Therefore, in order to continue our integrations further, we must 
consider a larger patch, such that the curve of 5 vs. m turns over. 
Our simulations do not show how much larger the patch must 
be, or indeed if any patch is adequate and A^-body simulations 
at later times must model a full annulus. Inconsistencies be- 
tween our A^-body integration and the Wetherill (1990) model 
rule it out to estimate the mass distribution. We note, how- 
ever, that more comphcated models {e.g., Kenyon and Bromley 
2004) may make a better match to our calculations, but the de- 
velopment and implementation of such models was beyond the 
scope of this investigation. We encourage future statistical re- 
searchers to use our results to verify their collision kernels. 
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Fig. 21 . — Evolution of all RMS velocities. 

In Fig. 21 (see also Figs. 7 and 19) we examine the differ- 
ent evolutions of the RMS velocity dispersion for all our sim- 
ulations. In each case, viscous stirring increases the velocity 
dispersion. Moreover, all have vrms « 2 m s~^ at ?i/2- Note 
that at fi/2 the average mass is 2m 1, which corresponds to an 
escape speed of 1 .7 m s~' . This equivalence suggests that, at 
least early on, the velocity dispersion grows at approximately 
the same rate as the escape speed of the typical mass particle. 

The perfect accretion model has produced several results that 
may be spurious. The sum of all these issues, in our imple- 
mentation, is that orbital angular momentum is too easily trans- 
ferred into rotational angular momentum, and produces spin 
rates that are too high (see Figs. 1 1 and 20). Therefore the spin 
distributions presented here should be not be regarded as phys- 
ically realistic, and the mass distributions of our simulations 
should be considered upper bounds. 

Assuming the initial swarm of planetesimals is composed of 
gravitational aggregates, then we need larger simulations with 
more a realistic coUisional model (see Leinhardt and Richard- 
son 2005) in order to determine on the true nature of the post- 
colhsion particles. Presumably angular momentum is lost by 
shedding rubble from the surface. Three possibilities await this 
freed rubble: collision with other planetesimals, orbital migra- 
tion via gas drag into the central star, or recollapse into 1 km 
planetesimals via collapse (Goldreich et al. 2004). Given the 
number density of planetesimals at this stage of growth, the for- 
mer seems the most likely. This linaitation of our model demon- 
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strates the need to perform similar simulations with a more re- 
alistic collisional model (e.g., Leinhardt and Richardson 2005; 
Leinhardt e/ a/. 2009). 

Our analytic fits to the mass distributions show that expo- 
nential fits do not match the data (see Table 3). However the 
power law does provide a reasonable fit for most of our mod- 
els. Removing the k = 1 bin from our fits results in a significant 
decrease in the unreduced values, as shown in Table 3. We 
conclude that runaway growth has not begun in our simulations. 

All of our simulations show that growth from 1 km planetesi- 
mals can proceed quickly. In fact, it proceeds so rapidly that our 
model breaks down in just a few hundred orbits (see §4.4.2, Fig. 
15). Therefore the only way to realistically proceed beyond ti/2 
is to enlarge our patches. As we expand our patches, the num- 
ber of particles increases. Given that Simulation Li required 
~ 30,000 node hours to complete on the Columbia Supercom- 
puter at NASA Ames, and computation time scales as Mog A'^, 
we may not be able to expand our patches such that they are 
both statistically accurate and computationally tractable. Our 
results therefore suggest that the patch model may be inade- 
quate to model later stages of terrestrial planet formation. 

7. CONCLUSIONS 

We have performed the first A'-body simulations of growth 
from 1 km planetesimals. The initial conditions of our runs 
were chosen to be similar to those believed to have existed in 
our protosolar disk; substantially different parameters may be 
appropriate for the initial stages of growth in extreme exoplan- 
etary systems (Lissauer & Slartibartfast 2008). These simu- 
lations required hundreds of thousands of node hours of su- 
percomputer time, using an advanced A'-body code designed 
specifically to examine systems with large A^. Although nu- 
merous shortcuts and approximations were incorporated in our 
model, we beheve that our results provide insights into plan- 
etesimal growth and lay a foundation for future investigations. 

Planetesimal growth from a uniform swarm of 1 km-sized 
planetesimals proceeds in a stochastic fashion. Our results have 
confirmed some of the trends seen in the semi-analytical re- 
search (Greenberg et al. 1978; Wetherill 1990; WeidenschiUing 



et al. 1997) into the growth of 1 km planetesimals. Until more 
realistic models of fragmentation can be implemented in A^- 
body code, statistical methods are the only feasible approach 
to address fragmentation. 

Some of the assumptions of our model broke down relatively 
quickly, demonstrating the limits of the patch approximation in 
modeling planetary accretion. Nonetheless, our results suggest 
new directions of research for this epoch of planet formation. In 
particular, a more realistic collisional model (one in which ad- 
ditional small particles carry away excess angular momentum, 
i.e., fragmentation) seems most important. Such a model may 
suppress growth (as well as ehminating unphysical spins), and 
hence planetesimals would not grow so quickly. However, this 
approach is considerably more complex and numerically inten- 
sive than those presented here. Until these modifications can 
be made, our results represent the most accurate model of 1 km 
planetesimal growth available. 
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APPENDIX A: LIST OF SYMBOLS AND ABBREVIATIONS 
A = arbitrary constant 

Aij = coUisional probability coefficient in coagulation theory 

a = semi-major axis 

B = arbitrary constant 

b = power law fit parameter 

b' = power law fit parameter with m = mi bin excluded 

C = arbitrary constant 

c = exponential fit parameter 

c' = exponential fit parameter with m = m\ bin excluded 
D = arbitrary constant 
e = eccentricity 

srms = root mean square eccentricity of particles in a patch 

Fg = gravitational focusing factor 

/ = fraction of particles relative to initial number 

G = Newton's gravitational constant 

i = incUnation 

iRMS = root mean square inclination of particles in a patch 

j = counter in coagulation equations 

k = ratio of particle mass to mass of 1 km planetesimal 

kR = mass of a runaway particle relative to a 1 km planetesimal 

Li = Largest baseline simulation 

/ = counter in coagulation equations 

Mo.5 = simulation with the initial velocity dispersion magnitude 

set to half that of Simulation Li 

Ml = medium- sized baseline simulation 

M2 = simulation with the initial velocity dispersion magnitude 
set to twice that of Simulation Li 
Mann = mass in an annulus of the protoplanetary disk 
Mpatch = total mass inside a patch 

M0 = mass of the Sun 

= mass of the Earth 
m = mass 

mi = mass of 1 km planetesimal 

merit = particle mass at which, in one orbit, it colhdes with an 
equal mass of gas 

nimax = largest mass in a simulation 

ffhutx-i = second largest mass in a simulation 

fnmax-A = fifth largest mass in a simulation 

mmax-9 = tenth largest mass in a simulation 

rupi = mass of planetesimal 

< m >= average mass of planetesimals 

= number of bodies in a simulation 
A'o = initial number of planetesimals in a simulation 
Npatch = number of bodies in a patch 
Nk = number of bodies in a mass bin 
rik = number density of particles in mass bin k 
P = orbital period 
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Ppeak = peak of the spin period distribution 

Si = Simulation with same initial properties as Li, but only 

1/16 the size 

PlAB = Particle-in-a-box 

R = planetesimal radius 

Rmax = radius of largest planetesimal 

r = heliocentric radius 

rpatch = hehocentric radius of the center of a patch 
f = time 

S = stirring power of largest mass relative to that of aU other 
bodies 

S 1 = smallest-sized baseline simulation 

f = time 

fi/2 = the time required to reduce the total number of particle at 
distance r by 2 

those = longest timestep in a simulation 

tcross = crossing time for two planetesimals 

tntin = minimum timestep in a simulation 

u = center-of-mass speed of a patch in the x-direction 

V = volume 

V = velocity 

Vese = escape speed of a planetesimal 
vrms = root mean square speed of a patch 
Vx = speed in x-direction 
Vy = speed in y-direction 
= speed in z-direction 
W = size scale of a patch 

w = center-of-mass speed of a patch in the y-direction 
X = Cartesian coordinate that mimics heliocentric distance 
Xg = x position of guiding center of a planetesimal's epicycle 
y = Cartesian coordinate that mimics azimuthal position 
yg = y position of guiding center of a planetesimal's epicycle 
Z = height above/below midplane 
Zq = scale height of planetesimal disk 
/? = radial excursions of a planetesimal due to eccentricity 
Pmax = largest radial excursion of a particle in a patch 
/3max-\ = second largest radial excursion of a particle 
Pmax-i = third largest radial excursion of a particle 
/3max-3 = fourth largest radial excursion of a particle 
Pmax-A = fifth largest radial excursion of a particle 
( = number of rungs in a simulation 
T] = scale factor to determine timesteps 

8 = maximum apparent size of a cell for PKDGRAV to only 
use the hexadecapole moment 

9 = azimuthal position of a planetesimal in hehocentric coordi- 
nates 

ui = coUisional probabihty coefficient in constant solution of 
coagulation equation 
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V2 = collisional probability coefficient in linear solution of co- 
agulation equation 

P3 = collisional probability coefficient in product solution of 

coagulation equation 

p = volume mass density 

po = volume mass density at midplane 

ppi = mass density of a planetesimal 

S = surface density 

So = coefficient that scales surface density of planetesimal disk 

a = physical cross-section of a planetesimal 

api = gravitationally enhanced cross-section of a planetesimal 

< (^pi >= average gravitationally enhanced cross-section of 

planetesimals 

T = mean free time between planetesimal physical collisions 
(p = gravitational potential 

xl = unreduced value for power law fit to mass distribution 
xl' = unreduced x^ value for power law fit to mass distribution 
with m = m\ mass bin excluded 

Xf = unreduced x^ value for exponential fit to mass distribution 
}Cc' = unreduced value for exponential fit to mass distribu- 
tion with m = mi mass bin excluded 
^patch = Keplerian orbital frequency of a patch 

= vertical frequency due to Keplerian motion and the mass 
of the disk 



APPENDIX B: THE COAGULATION EQUATION 

Here we summarize the basics of the coagulation equation as 
presented by Wetherill (1990). The discrete form of the coagu- 
lation equation is 



dNk . 



ijniiij 



l+j=k 



1=1 



(Bl) 



where the generic indices j, k and / are just the ratio m/mi. In 
Eq. (B 1), Nk is the number of particles in bin k, Aji is the prob- 
ability of collision, nj is the total number density of particles 
of mass j, m is the mass, and mi is the mass of a 1 km plan- 
etesimal. The first term is the probabiUty of all combinations of 
particles of mass I and j that sum to equal the current mass bin 
k. The mass bin is the quantum of the mass spectrum. The fac- 
tor of 1/2 prevents the summation from counting all collisions 
twice (when / = j). The second term is the loss of particles from 
mass bin k to larger mass bins. Note that, despite the discrete 
nature of Eq. (Bl), for f > 0, it can predict a fractional number 
of particles in each bin. 

The collisional probabihty coefficient, Aij, is a function of 
the relative velocity of particle / to j, their masses, the number 
density of each bin, and the volume being considered. There- 
fore 

Aij=Aij(Vrei,mi,mj,ni,nj), (B2) 

where «/ (rij) is the number density of particles with mass / (j). 
Three forms of this function have been examined. The simplest 
solution is 

Aij = uu (B3) 
a constant. For linear dependence we assume 



A,j = 1^2(1 + j), 



(B4) 



a constant times the sum of the masses, and the dependence 
on velocities and densities has been subsumed into U2. These 



two possible solutions both fall under the category of orderly 
growth. A third solution, which is proportional to the product 
of the masses, assumes 



Aij = uslj. 



(B5) 



At any given orbit, Wetherill also gives equations for the colU- 
sional probabiUties as a function of time: 



and 



ui = 



^2 = 



2(1-/) 
Noft ' 

log/ 
Not ' 

2(1-/) 

Not 



(B6) 



(B7) 



(B8) 



The solutions to the constant and sum forms are 

Nk = Nof\\-ff-\ 

and 

N,=Nj'^f{l-fr'e-'('-f\ 
k\ 

respectively, where A'o is the initial number of particles, and / is 
just the fraction of the number of particles remaining at time t. 



(B9) 



(BIO) 



/ = 



No No 

The product solution to the coagulation equation is 

N,=No^^(l-fr'e-"'^'-f>, 



(BID 



(B12) 



and yields runaway growth. This however leads to the natu- 
ral problem that a runaway particle is a special particle, and it 
should not be treated as typical. This marks the breakdown in 
the PIAB model. These solutions to the coagulation equation 
are used in §4. In the text we refer to Eq. (B9) as constant coag- 
ulation, Eq. (BIO) as linear coagulation, Eq. (B12) as product 
coagulation. 

APPENDIX C: THE BASELINE MODEL TO 2000 ORBITS 

In the spirit of Icarus' flight to the Sun, we present results for 
the baseline simulations from tif2 to 2000 orbits here. As shown 
in §4.4.2, after ?i/2 larger mass bodies may significantly alter 
the velocity dispersion in the patches. However, the locations 
of such particles relative to the patch are unknown. The synodic 
period across the radial width of the patch is about 2400 orbits. 
We may therefore presume that by 2000 orbits, a large mass has 
entered the patch and significantly altered the dynamics, but the 
time and magnitude of the changes are unknown. Although the 
results in this appendix suffer from significant inconsistencies, 
we nonetheless present them here, as they represent the only 
A^-body simulation of growth from 1 km planetesimals to date. 
Table 4 lists some of the properties of the baseline models at 
2000 orbits. Results in this appendix should not be regarded as 
physically realistic simulations of planetesimal growth beyond 
hi2'. 

Figure 22 presents the mass distribution of particles in a for- 
mat similar to Fig. 2. At 2000 orbits, three particles in Li and 
one in Mi have reached masses larger than 5000 mi. By com- 
parison, the fourth largest planetesimal in Li is one-sixth as 
massive, 814 m\, see Fig. 24. 
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Fig. 22. — Mass spectrum at various times during the evolution of 
the baseline patches in log-log format. Note that the Mi simulation 
has been offset by 0.1 and Si by 0.2. 

Figure 23 presents the mass functions at 2000 orbits, along 
with the associated analytic fits to the data. At 2000 orbits the 
power law has become a better fit to the data than at t\ /2- Table 
5 lists the fit parameters (and corresponding measures of good- 
ness of fit) at 2000 orbits. The values of the key parameters 
b and b' all cluster near 1.9. At 2000 orbits the exclusion of 
the k=\ bin does not cause dramatic changes in the the unre- 
duced values, which is not surprising since only 10% of the 
particles remain 3Xk= 1. 




5 10 
logs (k) = logj (m/m,) 

Fig. 23. — The mass function of Simulations Li (top), Mi (middle) 
and Si (bottom) and the associated power-law fits to the data; see Eq. 
(21) and Table 5. 

Figure 24 shows the evolution of some of the largest particles 
in the Li patch. At 713 orbits, the two largest particles in the 
patch merge to form a 1610 wi mass object. This merger has 



important consequences for our assumptions about the statisti- 
cal accuracy of our patch, as shown below. 
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Fig. 24. — Evolution of some of the most massive particles in the Li 
run. Note the merger of two ~ 800mi objects at 713 orbits. The three 
largest particle end up with masses in excess of 5000 mi . 

In Fig. 25 we show the evolution of vrms compared to the 
escape speed of the largest and typical particle in Simulation 
Li. After 500 orbits, vrms appears to grow linearly, while v™' 
begins to level off (except for the major merger event at 716 
orbits). These qualitatively different growth rates suggest that 
the larger particles are beginning to significantly heat the patch. 
Note that even this dynamical heating is probably an underes- 
timate, as particles outside the patch with larger mass should 
have sheared into this patch by 700 orbits. 
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Fig. 25. — Evolution of the escape speed of the largest particle 
(dashed line), RMS velocity (solid line), and escape speed of the aver- 
age mass particle (dotted line) for Simulation Li. 

We compare RMS velocity and mass growth for the three 
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baseline runs in Fig. 26. The RMS velocities remain close to 
each other up to ^ fi/2 (top panel), but subsequently diverge. 
This divergence corresponds with the masses of the largest par- 
ticles in Li and Mi reaching 100 mi. Note that Simulation 
Ml has the highest velocity dispersion from 615 to 1939 orbits. 
This feature occurs because the largest particle in Mi is nearly 
as large as that in Li (middle panel), but since the Mi patch is 
smaller than Li , the largest particle in Mi contains a larger frac- 
tion of the patch mass (bottom panel), and is therefore a more 
effective stirrer (see below). Although the largest particle in 
Simulation Si contains approximately the same fraction of the 
total mass as the particles in the other patches, its actual mass is 
considerably smaller, and, hence, the velocity dispersion in Si 
remains lower than the others. 
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Fig. 26. — Top: Evolution of the RMS velocities for the three baseline 
models. The values for Simulations Li and Mi stay relatively close to 
each other, but the S i value remains lower. Middle: Comparison of the 
growths of the largest particles in the baseline simulations. The largest 
particles in Li and Mi remain within about a factor of 2 each other, 
while the largest particle in Si lags by about an order of magnitude 
after 1000 orbits. Bottom: Fraction of the patch mass absorbed into 
the largest particle. At all times, the fractions remain within a factor of 
a few of each other. 



In Fig. 27 we see the evolution of the three gravitational fo- 
cusing factors in the baseline simulations. All appear to grow, 
reach a maximum, and then decrease. The peak is at about 300 
orbits for Si and 600 for Mi. The largest run Li is double- 
peaked, at 400 orbits and 700 orbits. The turnover at 400 orbits 
in Li may be a statistical fluke that is dramatically corrected at 
716 orbits, or it may be that the turnover at 400 orbits is real, 
and that the event at 716 orbits is anomalous. The curves in 
Fig. 24 suggest the former. Att = 200 the largest particle in Li 
begins to grow significantly faster than the second largest, even 
though the two have approximately the same mass. Then at 
t = 400, the difference between the two begins to decrease sud- 
denly. The two are almost identical at f = 700 when two large 
particle merge at 716 orbits. At this point, the largest particle 
once again becomes significantly larger. 
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Fig. 27. — Evolution of Fg as a function of time for the three baseline 
simulations. 

In Fig. 28 we show the evolution of crms/ siniuMS- After ini- 
tially dropping from 2.35 (see Fig. 9), the value in Li remains 
close to 2.1 for the duration of the simulation, suggesting the 
initial drop is a transient effect. However, in Mi and Si, the 
evolution is steady from 200 to 600 or 900 orbits, respectively, 
and then appears to drop monotonically at later times. The over- 
all drop represents about a 10% change. 
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Fig. 28. — The ratio of crms to sinzaw.s as a function of time. The 
ratio in Li remains slightly over 2 for the duration of the simulation, 
but the ratio for Mi and Si trends down. 

In Fig. 29 we examine dynamical friction in the baseUne runs 
at 2000 orbits. We saw in Fig. 10 that there was a general trend 
of decreasing velocity with increasing mass, although the slope 
of this trend is so shallow that kinetic energies trend higher with 
increasing mass. The same general pattern is seen at 2000 or- 
bits for large masses (m >100mi), but velocity is independent 
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of mass for smaller masses. Note as well that the mean values 
for the small-mass particles are larger at 2000 orbits than at ?i/2- 
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Fig. 29. — The values of eRMs and Irms as a function of mass at 2000 
orbits for the three baseline simulations. Filled squares are eccentric- 
ity, open triangles are inclination. For reference, the dashed line repre- 
sents equipartition of energy in e, dotted in i, normalized to values at 
k= 10. 

Next we examine the validity of the Li patch from - 
2000 orbits. First we plot the evolution of u and w in Fig. 30. 
Both values remain small until the major merger event at 713 
orbits. At that time, the values grow significantly, finally reach- 
ing M = -7 cm s-\ about two orders of magnitude larger than its 
value at f 1 /2 ■ This velocity means that the center-of-mass veloc- 
ity is roughly 1/70 the total shear across the patch, and that cer- 
tainly by the end of the simulation our assumptions have broken 
down. 
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Fig. 30. — Evolution of u and w in Simulation Li. The values grow 
quickly after the merger at 713 orbits. 



Next we look at how effective different mass bins are at stir- 
ring the patches at different times in Fig. 31. We saw in Fig. 14 
that by t\j2 the distribution of m^Nk was flat, and in Fig. 31 we 
see that the slope is positive at all times after ?i/2 for all simula- 
tions. These positive slopes indicate that the patch is not large 
enough since very large particles not in the patch can substan- 
tially affect the velocity distribution in the patch (see § 4.4.2). 
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Fig. 3 1 . — Distribution of stirring power as a function of mass at four 
separate times in all baseline simulations. By the end of the simula- 
tion, each of the three largest particles in Li (solid line) are an order 
of magnitude more effective at stirring the patch than all the other par- 
ticles combined. Note that the Mi Simulation data (dashed line) are 
offset by 0.1 and the Si data (dotted Une) by 0.2. 



Finally we examine the value of S in Fig. 32. The value in 
Li stays below 0.15 until 713 orbits, indicating the largest-mass 
particle is not dominating the stirring in the patch. However, the 
merger at 713 orbits creates a particle whose stirring is about 
equal to the stirring of the sum of all other particles in the patch. 
Therefore at this point, the assumptions of the patch framework 
break down, and we can not expect the simulation to be pro- 
viding reliable results. As other large particles appear in Li, 5 
slowly drops. In Mi, the value nearly reaches unity at 700 or- 
bits, and then grows quickly to a final value of 17.5. In Si, 5 
remains below unity for the duration of the simulation, but note 
the sudden jump at 1600 orbits. 
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Fig. 32.- 
models. 



■ Evolution of S as a function of time for the three baseline 



The results presented in this appendix reveal the numerous 
ways in which our patch model breaks down after fi/2. The as- 
sumptions of small center-of-mass motion, no dominant mass 
inside the patch, and no outside perturbers have all failed in the 
three baseline models. 



Table 1. Initial Conditions of the Patch Simulations 



ID 


N 


Vrms (m S'O 


Z() (km) 


po (10-' g/cm^) 


Li 


106,130 


1.29 


475 


3.1 


Ml 


26,532 


1.29 


475 


3.1 


Mo.5 


26,532 


0.65 


237 


6.3 


M2 


26,532 


2.6 


950 


1.6 


Si 


6633 


1.29 


475 


3.1 


Table 2. 


Results of 1 km Planetesimal Growth at 0.4 AU at fi/2 


ID hi 


2 (orbits) 


m,„ax (mi) 


(m S" 


) Fg Pp,ak (hr) 


Li 


354 


276 


1.93 


20.2 1.05 


Ml 


349 


142 


1.89 


13.7 1.05 


Si 


361 


73 


1.95 


9.2 1.05 


Mo.5 


253 


202 


1.86 


17.8 1.65 


M2 


542 


47 


2.41 


4.8 0.65 



Table 3. Fit Parameters of the Mass Distributions at f 1/2 



ID 


b 


xl 


b' 


xi' 


c 




c' 


xl' 


Li 


2.39 


1693 


2.62 


376 


1.3 


6224 


1.9 


1690 


Ml 


2.39 


894 


2.64 


168 


1.2 


1400 


1.9 


853 


Si 


2.38 


492 


2.65 


107 


1.3 


1073 


1.9 


308 


Mo.5 


2.4 


307 


2.59 


74 


1.2 


1723 


2 


466 


M2 


2.38 


489 


2.61 


154 


1.3 


1440 


1.7 


189 



Table 4. Results of 1 km Planetesimal Growth at 0.4 AU after 2000 orbits 



ID 


m„uvc im) 


Vrms (m S"') 




Ppeak (hr) 


Li 


5879 


7.94 


9.69 


0.85 


Ml 


4117 


7.89 


7.95 


0.75 


Si 


113, 


6.00 


4.95 


0.85 



Table 5. Fit Parameters of the Mass Distributions at 2000 Orbits 



ID 


b 


Xb 


b' 


Xb' 


c 


xi 


c' 


xl' 


Li 


1.88 


486 


1.93 


253 


8.0 


11562 


10.2 


4511 


Ml 


1.89 


141 


1.96 


71 


7.4 


2783 


9.8 


1029 


Si 


1.85 


31 


1.91 


20 


7.9 


560 


10.0 


206 



